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Abstract. 

This paper gives a detailed pedagogic presentation of the central con- 
cepts underlying a new algorithm for the numerical solution of Einstein's 
equations for gravitation. This approach incorporates the best features of 
the two leading approaches to computational gravitation, carving up space- 
time via Cauchy hypersurfaces within a central worldtube, and using char- 
acteristic hypersurfaces in its exterior to connect this region with null infin- 
ity and study gravitational radiation. It has worked well in simplified test 
problems, and is currently being used to build computer codes to simulate 
black hole collisions in 3-D. 



1. Preamble 

Throughout his career in theoretical gravitational physics, Vishu has been 
interested in questions of black holes and gravitational radiation. The sub- 
ject of his useful early study of radiation from binary systems [1] is now at 



2 



NIGEL T. BISHOP ET AL. 



the forefront of research as the principal target for the first exciting experi- 
mental measurements of the LIGO project. His fundamental studies of the 
properties of black hole excitations and radiation [2] used the formidable 
technology of the mid 1960's, i.e. Regge- Wheeler perturbation theory, to 
analytically extract the crucial physical result that black holes were stable. 
In his continuing studies of the interaction of black holes with gravitational 
radiation [3], he first demonstrated the phenomenon that was later to be 
called normal mode excitations of black holes, and noted that the frequency 
of the emitted radiation carried with it key information which could be used 
to determine the mass of the invisible black hole. These issues are still at 
the heart of current research three decades later. 

Today, the technology to study these questions has become even more 
formidable, requiring large groups of researchers to work hard at devel- 
oping computer codes to run on the massively parallel supercomputers of 
today and tomorrow. This has caused gravitational theory to enter into the 
realm of "big science" already familiar to experimental physics, with large, 
geographically distributed collaborations of scientists engaged on work on 
expensive, remote, central facilities. The goal of this modern work is to un- 
derstand the full details of black hole collisions, the fundamental two-body 
problem for this field. Recent developments in this area are very encourag- 
ing, but it may well take another three decades until all the riches of this 
subject are mined. 

In this paper, we will present all the gory details of how the best current 
methods in computational gravitation can be forged into a single tool to 
attack this crucial problem, one which is currently beyond our grasp, but 
perhaps not out of our reach. 

2. Introduction 

Although Einstein's field equations for gravitation have been known for the 
past 80 years, their complexity has frustrated attempts to extract the deep 
intellectual content hidden beneath intractable mathematics. The only tool 
with potential for the study of the general dynamics of time-dependent, 
strongly nonlinear gravitational fields appears to be computer simulation. 
Over the past two decades, two alternate approaches to formulating the 
specification and evolution of initial data for complex physical problems 
have emerged. The Cauchy (also known as the ADM or "3 + 1") approach 
foliates spacetime with spacelike hyper surfaces. Alternatively, the charac- 
teristic approach uses a foliation of null hyper surf aces. 

Each scheme has its own different and complementary strengths and 
weaknesses. Cauchy evolution is more highly developed and has demon- 
strated good ability to handle relativistic matter and strong fields. However, 
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it is limited to use in a finite region of spacetime, and so it introduces an 
outer boundary where an artificial boundary condition must be specified. 
Characteristic evolution allows the compactification of the entire spacetime, 
and the incorporation of future null infinity within a finite computational 
grid. However, in turn, it suffers from complications due to gravitational 
fields causing focusing of the light rays. The resultant caustics of the null 
cones lead to coordinate singularities. At present, the unification of both of 
these methods [4] appears to offer the best chance for attacking the funda- 
mental two-body problem of modern theoretical gravitation: the collision 
of two black holes. 

The basic methodology of the new computational approach called Cauchy- 
characteristic matching (CCM), utilizes Cauchy evolution within a pre- 
scribed world-tube, but replaces the need for an outer boundary condi- 
tion by matching onto a characteristic evolution in the exterior to this 
world-tube, reaching all the way out to future null infinity. The advantages 
of this approach are: (1) Accurate waveform and polarization properties 
can be computed at null infinity; (2) Elimination of the unphysical out- 
going radiation condition as an outer boundary condition on the Cauchy 
problem, and with it all accompanying contamination from spurious back- 
reflections, consequently helping to clarify the Cauchy initial value problem. 
Instead, the matching approach incorporates exactly all physical backscat- 
tering from true nonlinearities; (3) Production of a global solution for the 
spacetime; (4) Computational efficiency in terms of both the grid domain 
and algorithm. A detailed assessment of these advantages is given in Sec. 3. 

The main modules of the matching algorithm are: 

— The outer boundary module which sets the grid structures. 

— The extraction module whose input is Cauchy grid data in the neigh- 
borhood of the world-tube and whose output is the inner boundary 
data for the exterior characteristic evolution. 

— The injection module which completes the interface by using the ex- 
terior characteristic evolution to supply the outer Cauchy boundary 
condition, so that no artificial boundary condition is necessary. 

Details of the Cauchy and characteristic codes have been presented else- 
where. In this paper, we present only those features necessary to discuss 
the matching problem. 

3. Advantages of Cauchy-characteristic matching (CCM) 

There are a number of places where errors can arise in a pure Cauchy 
computation. The key advantage of CCM is that there is tight control over 
the errors, which leads to computational efficiency in the following sense. 
For a given target error e, what is the amount of computation required for 
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CCM (denoted by Aqcm) compared to that required for a pure Cauchy 
calculation (denoted by Awe)? It will be shown that Accm/Awe — * O 
as e — > O, so that in the limit of high accuracy CCM is by far the most 
efficient method. 

In CCM a "3 + 1" interior Cauchy evolution is matched to an exterior 
characteristic evolution at a world-tube of constant radius R. The impor- 
tant point is that the characteristic evolution can be rigorously compacti- 
fied, so that the whole spacetime to future null infinity may be represented 
on a finite grid. From a numerical point of view this means that the only 
error made in a calculation of the gravitational radiation at infinity is that 
due to the finite discretization h; for second-order algorithms this error is 
0{h 2 ). The value of the matching radius R is important, and it will turn 
out that for efficiency it should be as small as possible. The difficulty is that 
if R is too small then caustics may form. Note however that the smallest 
value of R that avoids caustics is determined by the physics of the problem, 
and is not affected by either the discretization h or the numerical method. 

On the other hand, the standard approach is to make an estimate of 
the gravitational radiation solely from the data calculated in a pure Cauchy 
evolution. The simplest method would be to use the raw data, but that ap- 
proach is too crude because it mixes gauge effects with the physics. Thus a 
substantial amount of work has gone into methods to factor out the gauge 
effects and to produce an estimate of the gravitational field at null infinity 
from its behavior within the domain of the Cauchy computation [5, 6, 7]. 
We will call this method waveform extraction, or WE. The computation is 
performed in a domain D, whose spatial cross-section is finite and is nor- 
mally spherical or cubic. Waveform extraction is computed on a world-tube 
r, which is strictly in the interior of D, and which has a spatial cross-section 
that is spherical and of radius te- While WE is a substantial improvement 
on the crude approach, it has limitations. Firstly, it disregards the effect, 
between V and null infinity, of the nonlinear terms in the Einstein equations; 
the resulting error will be estimated below. Secondly, there is an error, that 
appears as spurious wave reflections, due to the inexact boundary condi- 
tion that has to be imposed at dD. However, we do not estimate this error 
because it is difficult to do so for the general case; and also because it is in 
principle possible to avoid it by using an exact artificial boundary condition 
(at a significant computational cost). 

The key difference between CCM and WE is in the treatment of the 
nonlinear terms between T and future null infinity. WE ignores these terms, 
and this is an inherent limitation of a perturbative method (even if it is 
possible to extend WE beyond linear order, there would necessarily be a cut- 
off at some finite order) . Thus our strategy for comparing the computational 
efficiency of CCM and WE will be to find the error introduced into WE from 
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ignoring the nonlinear terms; and then to find the amount of computation 
needed to control this error. 

3.1. ERROR ESTIMATE IN WE 

As discussed earlier, ignoring the nonlinear terms between T (at r = t\e) 
and null infinity introduces an error, which we estimate using characteristic 
methods. The Bondi-Sachs metric is 

ds 2 = -^j-r 2 h AB U A U B y u 2 

- 2e 2l3 dudr - 2r 2 h AB U B dudx A + r 2 h AB dx A dx B , (1) 

where A, B = 2,3 and h AB is a spherical metric that is completely de- 
scribed by one complex function J. The initial data required on a null cone 
u = constant is J, and the hypersurface equations R\ a = then form 
a hierarchy from which (3, U and W = (V — r)/r 2 are found [11]. The 
evolution equation R AB h AB = [11] is 

2(rJ), ur = Lj + Nj (2) 

where Lj represents the linear part and Nj the nonlinear part. 

The order of magnitude of various terms can be expressed in terms of 
a function c(u, x A ) (whose time derivative c u is the news function); note 
that c is not a small quantity. The expressions are 

J = 0( C -),f3 = O(^), U A = 0(£), W = O(^). (3) 

These estimates are obtained from the hypersurface equations, and assume 
that the background geometry is Minkowskian. Should this not be the case 
then constants of order unity would be added, and the effect of this would be 
to amend (2) by adding terms to Lj so that it represents wave propagation 
on a non-Minkowskian background. However, the order of magnitude of 
terms in Nj would not be affected. It is straight forward to confirm that 
Nj involves terms of order 

Oi^)- (4) 

WE estimates the news at future null infinity from data at r = r B , and 
could be made exact if Nj were zero. Thus the error introduced by ignoring 

Nj is 

roc £,2 ^2 

e(c«) = (c u )exact ~ (c u )we = / 0(—)dr = O(-). (5) 

Jr E r r E 
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This would be the unavoidable linearization error in WE were it imple- 
mented as an "exact" artificial boundary condition by using global tech- 
niques, such as the difference potential method, to eliminate back reflec- 
tion at the boundary [8]. However, this is computationally expensive [9] 
and has not even been attempted in general relativity. The performance of 
WE continues to improve but the additional error due to back reflection 
remains [10]. 

3.2. COMPUTATIONAL EFFICIENCY 

A numerical calculation of the emission of gravitational radiation using a 
CCM algorithm is expected to be second-order convergent, so that after a 
fixed time interval the error is: 



where h is the discretization length. On the other hand, the same calculation 
using WE must allow for the error found in (5), and therefore after the same 
fixed time interval there will be an error of: 



We now estimate the amount of computation required for a given desired 
accuracy. We make one important assumption: 
— The computation involved in matching, and in waveform extraction, 
is an order of magnitude smaller than the computation involved in 
evolution, and is ignored. 
For the sake of transparency we also make some simplifying assumptions; if 
not true there would be some extra constants of order unity in the formulas 
below, but the qualitative conclusions would not be affected. 

1. The amount of computation per grid point per time-step, a, is the 
same for the Cauchy and characteristic algorithms. 

2. The constants k±,k2 in the equations above are approximately equal 
and will be written as k. 

3. In CCM, the numbers of Cauchy and characteristic grid-points are the 
same; thus the total number of grid points per time-step is: 



4. In WE, the number of grid points in D is twice the number contained 
in T; thus the total number of grid points per time-step is 



£ = Q{h 2 ) ~ hh 2 , 



(6) 



e = 0(h 2 ,r E 2 ) ~ maxfah 2 , -|-). 



(7) 




(9) 
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It follows that the total amount of computation A (i.e., number of floating- 
point operations) required for the two methods is: 

87ri? 3 a 8irr E a 
Accm = ^j^, A WE = ^ r . (10) 

Thus R > te or R < r E determines which method requires the least amount 
of computation. Because of the assumptions (1) to (4) this criterion is not 
exact but only approximate. 

As stated earlier, in a given physical situation the minimum allowed 
value of R is determined by the physics. However, is determined by the 
target error (equation (7)); and there is also a minimum value determined 
by the condition that the nonlinearities must be sufficiently weak for a 
perturbative expansion to be possible. Thus, in a loose sense, the minimum 
value of r E is expected to be related to the minimum value of R. It follows 
that the computational efficiency of a CCM algorithm is never expected to 
be significantly worse than that of a WE algorithm. 

If high accuracy is required, the need for computational efficiency always 
favors CCM. More precisely, for a given desired error e, equations (6) and 
(7), and assumption (2), imply: 



h = yje/k, r E = ^k 3 /e. (11) 
Thus, substituting equation (11) into equation (10), 



8nR 3 ak 2 Svrafc 2 ^ 72 

A C CM = ^ ' A WE = 3£7/2 : 

so that 



(12) 



Accm R 3 e 3/2 n . 
— = — ^ ► as e -> 0. (13) 

Awe kl /2 

This is the crucial result, that the computational intensity of CCM relative 
to that of WE goes to zero as the desired error s goes to zero. 



4. Extraction 

We describe a procedure by which information about the 4-geometry on 
the neighborhood of a world-tube T, obtained during a 3 + 1 simulation 
of Einstein's equations, is used to extract boundary data appropriate for 
an exterior characteristic formulation. A numerical implementation of this 
characteristic formulation [12], is then used to propagate the gravitational 
signal to null infinity, where the radiation patterns are calculated. The 
process we describe here is non-perturbative, and to make it as portable as 
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possible, it assumes only that the 3 + 1 simulation can provide the 3- metric, 
the lapse and the shift by interpolation to a set of prescribed points. 

We start by describing the world-tube T in Sec. 4.1, we explain what 
information needs to be provided by the 3 + 1 simulation around the world- 
tube in Sec. 4.2, giving details of the transformation to null coordinates in 
Sec. 4.3, and the expression for the Bondi metric in Sec. 4.5. This provides 
the boundary data needed to start up the characteristic code in Sec. 4.6. The 
characteristic code which takes this boundary information and calculates 
the waveforms is described in detail in [12]. 

4.1. PARAMETRIZATION OF THE WORLD-TUBE 

The notation and formalism are based on Misner, Thorne and Wheeler [13]. 
Greek indices range from 1 to 4, Latin indices range from 1 to 3, and upper 
case Latin indices refer to coordinates on the sphere and range from 2 to 3. 

The intersections St of the world-tube with the space-like slices of the 
Cauchy foliation are topologically spherical, and they can be parametrized 
by labels y A , A = [2,3] on the sphere. The intersections themselves are 
labeled by the time coordinate of the Cauchy foliation, x 4 = t. Future 
oriented null cones emanating from the world-tube are parametrized by the 
labels on the sphere y A and an affine parameter A along the radial direction, 
with A = on the world-tube. With the identifications y 1 = A and retarded 
time y 4 = u = t, we define a null coordinate system y a = (y 1 ,y A ,y 4 ). 
We will later introduce a second null coordinate system y a = (y 1 , y A , y 4 ) , 
where y 1 = r, r being a surface area coordinate, and y A = y A , y 4 = y 4 . 

Following [14], we cover the unit sphere with two stereographic coor- 
dinate patches, centered around the North and South poles, respectively, 
where the stereographic coordinate is related to the usual spherical polar 
coordinates (9,(j)) by 

/ l-cosfl i4> / l + cosfl i4> 

^North = \ T~. v , Csouth = \ r. ' p . (14) 

y 1 + cos 6 y 1 — cos 

Let y A = (q,p), for A = [2,3], label the points of a stereographic coordi- 
nate patch. We adopt as these labels the real and imaginary part of the 
stereographic coordinate, i.e. for £ = y 2 + iy 3 = q + ip. 

For the computational implementation of the procedure we are describ- 
ing, we introduce a discrete representation of these coordinate patches 

y2 = -i + (i-3)A, yJ = -l + (j-3)A, (15) 

where the computational spatial grid indices i , j , k run from 1 to N and 
A = 2/(N — 5), N is the stereographic grid size and A the stereographic 
grid spacing. 
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On each patch, we introduce complex null vectors on the sphere q A = 
(P/2)(5 A + iS A ), where P = 1 + ££. The vectors q A , q A define a metric on 
the unit sphere 



1 , \ 4 

9 AS = ^ WA9B + 9A9.b) = -pf 



1 
1 



(16) 



with determinant det(gAB) = 16/P 4 . The co- vector ^4 satisfies the orthog- 
onality condition q~AQ = 2, and has components qA = (2/-P) (S\ + iS\). 

Given Cartesian coordinates x % = (x,y,z) on a space-like slice T, t , the 
intersection S t of the world-tube T with Ej is described parametrically by 
three functions of the y^jX*- -** = f l (y A )- In the following, we will fix the 
location of the world-tube by choosing these functions (in stereographic 
coordinates) as : 



f x (y A ) = 2R 



f y (y A ) = ± 2R (^i) 



f z (y A ) = ±r (tt§) 

where the positive (negative) sign corresponds to the north (south) patch. 
This defines a canonical spherical section St of radius R in Minkowski-space. 
Note that this provides a prescription for locating the world-tube which is 
time-independent. 

Given two 3-dimensional quantities v l and w l , we introduce their Eu- 
clidean inner product (y ■ w) = 5ijv' l w J . The stereographic coordinates 
of points on the surface of the 3-sphere with Cartesian coordinates x l = 
(x,y,z) can then be represented by y A {x l ) = (q,p), with 

</(**) = ^-, p{x i )=±J- (18) 
r ± z r ± z 

on the north (+) and south (— ) patches, where 

f 2 = {x-x). (19) 

Note that the stereographic coordinates are invariant under a change of 
scale, y A {x % ) = y A (cx % ). 

The Euclidean radius of the world-tube with Cartesian coordinates x^ ^ 
is given by: 

R? = { x W-xW). (20) 



10 



NIGEL T. BISHOP ET AL. 



This equation provides a particularly simple alternate definition of the choice 
of world-tube using Cartesian coordinates. This definition also holds for 
each instant of Cartesian time. 

4.2. A-D GEOMETRY AROUND THE WORLD-TUBE 

The A-d geometry around the world-tube is fully specified by the 4-d metric 
and its derivatives (alternatively, by the metric and the metric connection) . 
They determine the unit normal n a to the t = constant slices, and given the 
parametrization of the world-tube, the (outward pointing) normal s a to the 
world-tube. They also determine the generator of the outgoing null radial 
geodesies through the world-tube, which in turn completes the specification 
of the coordinate transformation x a — > y a in a neighborhood of the world- 
tube. 

In practice, the necessary information is not available at the discrete 
set of points on the world-tube specified by Eqs. (14) and (17) where £ij = 
y\ + iy^ as given in Eq. (15). However, the required variables are known 
on the points of the computational grid used in the simulation, a Cartesian 
grid (xi, y-j , z k ), from which we interpolate them to the world-tube points 
to second order accuracy. 

For a standard 3 + 1 simulation [15], the variables that we need to in- 
terpolate are the 3-d metric gij, the lapse a and the shift Their spatial 
derivatives are also interpolated. Their values at the world-tube points arc 
stored for a number of time levels, and the time derivatives of the 3-metric, 
lapse and shift at the world-tube are computed by finite-differencing be- 
tween these time levels. 

Using all these values we can compute the 4- metric g^, and its first 
derivatives g^a, using the following relations: 

9it = 9ijf3 J 

9tt = -a 2 + gu(3 l 

9it,n = 9ij,^ 3 + 9ijP^ 

g tt ^ = -loo,, ■ u,,,, ■ 2 9ij f>„ . (21) 

The unit normal n M to the hypersurface St is determined from the lapse 
and shift 

n" = ±(1,-/3*). (22) 

Let s a = (s 1 , 0) be the outward pointing unit normal to the section St of 
the world-tube at time t n . By construction, s* lies in the slice S t , and it is 
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known given the two vectors dy2 , d y 3 in S t , defined by the parametrization 
of the world-tube x l {y A ) 

«' = W P ' = W <23) 

These may be obtained analytically from Eq. (17), the equation for the 
world tube. Antisymmetrizing q l and p l , we obtain the spatial components 
of the normal 1-form cr 9 ; and its norm a 



°i = tijk(l J V k , cr = yJg ij ai(Tj (24) 

from which s l is obtained by raising <7j with the contravariant 3-metric g lJ 
on the slice T, t and dividing by the norm a, yielding 

8* = ^^. (25) 

The generators £ a of the outgoing null cone Ct through St are given on 
the world-tube by 

n a 4- s a 

a- gij p l sJ y ' 

which is normalized so that i a t a = — 1, where t a = an a + (3 a is the Cauchy 
evolution vector. 

The equations in this section show explicitly how to use the output data 
from the Cauchy simulation to completely reconstruct the full 4-geometry 
of the spacetime, as well as other important geometrical objects of interest, 
in the neighborhood of the world-tube T. This is all described in this sec- 
tion within the Cartesian coordinate system used by the 3 + 1 computation, 
and holds to the second-order accuracy assumed for this computation. In 
the next three sections, we will demonstrate how to use this information 
to redescribe the same geometry in another coordinate system, the Bondi 
coordinates needed for characteristic simulations. There is nothing in these 
sections that goes beyond the elementary concepts of defining coordinate 
systems and transforming tensors under a change of coordinates. Neverthe- 
less, it is quite instructive to see just how much work lies hidden beneath 
the clever notation used by theorists. 



4.3. COORDINATE TRANSFORMATION 

In this section we build the coordinate transformation between the 3 + 
1 Cartesian coordinates x a and the (null) afhne coordinates y a . We will 
need this in the neighborhood of the world-tube, not just at a point on 
the tube, in order to easily pass information back and forth between the 
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Cauchy and characteristic computer codes in the overlap region near the 
world-tube. In Sec. 4.4, we will transform the metric to affine coordinates. 
Finally, in Sec. 4.5, we will complete the transformation from affine to 
Bondi coordinates. 

The motivation for this indirect route to the Bondi coordinate frame 
deserves some comment. Both affine coordinates y a and Bondi coordinates 
y a utilize null hypersurfaces for foliations. Calculating these directly would 
require the numerical solution of a nonlinear partial differential equation 
(the eikonal equation). A much simpler, but physically equivalent, approach 
is to instead solve the null geodesic equation in Cartesian coordinates, in 
order to find the rays x^(X) generating the required null hypersurfaces. 
Secondly, we must introduce the intermediary of the affine radial coordinate 
A rather than the Bondi surface area coordinate r, since the latter is actually 
unknown until the angular coordinates are defined. As shown below, it is 
only after the null rays have been found that we are able to proceed with the 
orderly introduction of angular coordinates y A in the exterior of the world- 
tube. Finally, null geodesies can be analytically constructed trivially in 
terms of A, whereas even when we have defined the surface area coordinate 
r, solution of the null geodesies equation requires the numerical solution of 
an ordinary differential equation. With this rationale complete, we will now 
proceed to carry out the coordinate definitions and transform the metric. 

By inspection, x a (X), the solution to the geodesic equation relating x a 
to y a off the world-tube is: 

x a = x (0)a + ^ (0)a A + O(A 2 ). (27) 

This expression determines x a (X) to 0(A 2 ), given the coefficients 

x ( ® a = x c \ T and £^ a = x% ]r (28) 

that is, given the coordinates of the points and the generators of the null 
cone through the world-tube section St- For completeness, we repeat the 
remaining coordinate relations defined in Sec. 4.1. Along each outgoing 
null geodesic emerging from St, angular and time coordinates are defined 
by setting their values to be constant along the rays, and equal to their 
values on the world-tube 

y A = y A r and f = u = t. (29) 

Given the coordinate transformation x^ = x^(y a ), we obtain the metric 
in null affine coordinates r) a p by 



cte M dx y 

'a/3 = d ya Qyfi^v- 



(30) 
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The Jacobian of the coordinate transformation is viewed as a series expan- 
sion in the affine parameter A for each point on the world-tube. Further- 
more, we may omit the A derivatives of the x^, i.e. the are not needed 
because the radial coordinate A is an affine parameter of the null geodesies, 
hence the fjxp, components of the null metric are fixed: 

^?AA = VxA = °> Vxa = -1, (31) 

(This follows from the conditions s a n a = 0, £ a £ a = 0, s a s a = 1, n a n a = 
— 1 and t a £ a = —1). In other words, fj & ^ contains six independent metric 
functions. The relevant part of the coordinate transformation is then 



x%^^ = x^% + x^%X + 0(X% xM% = l^%, for a = (A,u). 

(32) 

Because the specification of the world-tube is time-independent, and the 
slices St of the world-tube are by construction at t = constant, only the 
angular derivatives of the x^ 1 for i = 1,2,3 survive, i.e. the O(A ) part of 
the Jacobian is given by the condition di/dui^ = 1, and by the relations 



which can be evaluated from the analytic expressions Eq. (17). 
To evaluate the 0(A) part of the Jacobian, we note that 

x M = i ^ X M = Z% ( 34 ) 
We will spend the remainder of this section evaluating these terms. To 
proceed, we see from Eq. (26) that this will require derivatives of n M and 
s l . For the simple case of geodesic slicing, a = 1, j3 l = 0, the derivatives 
of the normal n^ vanish, and the transformation is contained purely in the 
derivatives of s^. In general, the angular derivatives of n^ at A = constant 
can be computed in terms of spatial 3+1 derivatives of the lapse and shift, 
transformed with the 0(A°) Jacobian, i.e. 



The 3+1 derivatives are given by 



n 



n 



a 



(36) 
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The retarded time derivative da at A = constant is simply the 3 + 1 time 
derivative dt, therefore = n^t where 



n\ t = «t- (37) 

From Eq. (25), and since the G k are time-independent, the time deriva- 
tive of s 1 is given by 



s ,t - 9,t a 9 a 2 - ~9 9 9mn,t a ~ s a 

= -g tm 9mn,tS n -S^ (38) 

a 

where the time derivative of a can be calculated from 



o / 2\ kl km In m n 2 /on\ 

2aa tt = [a ) f = 9,t °k<?i = ~9 9 g m n,t<Jk<?i = -s s g mn ,to , (39) 
with the resulting expression 

s\ t = [-9 im + s % \s m ^g mn , t s n . (40) 
Similarly, it follows from Eq. (25) that 



J _ Jk T j °k ik a k,A ik akU ,A 

,A U ,A (J (J (J* 



g m g km g mn , ^ p- + g lk ^ - ^ (41) 



> A g a a 



where the G k A are obtained from the analytic expressions Eq. (17), and 
g ^ from 

2gg a = (a 2 ) j = (g kl G k a l ) ^ = g^aw + 2 / W M 

km In j , o kl 

= ~g g 9mn,jXy k Gl + 2 g G t G kA 

= - s m s n g mn , j x :i A G 2 + 2s k GG k A (42) 

Collecting Eqs. (41) and (42), we arrive at the angular derivatives of 
the normal to the world-tube 
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jn „n „ n k K ,s± 



s J[ — 9 $ 9mn,j % 9 g "1 ~ $ \^2 S ^ 9mn,j% S 

= (gin _ s . s n) °J±A + (-gin + ^ A ^.jj (43) 



4.4. NULL METRIC f\ aSi 

Given the 4-metric and its Cartesian derivatives at the world-tube, we can 
calculate its derivative with respect to the affine parameter A according to 

Having obtained the relevant parts of the coordinate transformation 
y a -» x a , Eqns. (34), (35), (36), (37), (40), and (43), and given the metric 
and its A derivative as per Eq. (44) , we can expand the null metric as follows 

^ = €j+%/3,AA + 0(A 2 ), (45) 
where the coefficients are given by 



Vm = 9tt\r 

~(0) i 

\A = X ,A3it\T 



'AB ,A ,B 

and, for the A derivative 



€b = *'i^m- ( 46 ) 



Vuu,\ 
VuA,\ 

Vab,\ 



9tt,x + 2£%g, t , +0(A) 



( t~ 

A V ,« 



u 9kfi + 9kt,x) + ^9kt + ^x9tt 



|r 



+ 0(A), 



x k A x l ^g k i,x + 



^Ab + ^a 



9 id 



+ 0(A). (47) 



The remaining components are fixed by Eq. (31). 

It is worthwhile to discuss a subtlety in the rationale underlying the 
computational strategy used here. The purpose of carrying out an expan- 
sion in A is to enable us to give the metric variables in a small region off the 
world-tube, i.e. at points of the grid used to discretize the null equations. 
This method is an alternative to the more obvious and straightforward 
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strategy of interpolation to determine metric values at needed points on 
the null computational grid. At first glance, it might appear a more cum- 
bersome way to proceed, due to the necessity of analytic calculation of the 
derivative of the metric on the world tube to determine coefficients needed 
for the metric expansion in powers of A. On the contrary, however, this trick 
allows a major simplification in computer implementation. It allows us to 
use a null code which does not require any special implementation at the 
irregular boundary defined by the world-tube, and it automatically ensures 
continuity of the metric and the extrinsic curvature at the world-tube. The 
obvious alternative of a brute-force interpolation approach would require a 
full 4-dimensional evaluation with quite complicated logic, since it would lie 
at the edge of both the Cauchy and null computational grids. Instead, the 
A expansion reduces the complexity by one dimension, allowing for much 
easier numerical implementation. 

We also need to compute the contravariant null metric, r/ a/3 , which we 
similarly consider as an expansion in powers of A, 



It follows from Eq. (31) that the following components of the contravari- 
ant null metric in the y a coordinates are fixed 



(48) 



with coefficients given by 



~(0)/ia~(0) _ <rjj. ~hv _ _~fia~f3v ~ _ 
'/ 'lav — °ii ' 'l,X — 'I 'I 'la/3 



(49) 



7? A « = -1 fj uA = fj UU = 0, 

therefore the contravariant null metric can be computed by 



(50) 





-Vuu + V VAu 



(51) 



and similarly for its 



A derivative 





= -Vuu,x + 2 r A %A,x - V V Vab,x 



(52) 
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4.5. METRIC IN BONDI COORDINATES 

The surface area coordinate r(u, A, y A ) is defined by 

r ~\det{q AB )) ' (53) 

where, for our choice of stereographic coordinates £ = q + ip, we use y A = 
(q,p) and det(qAB) = 16/(1 + q 2 + p 2 ) 4 • To carry out the coordinate 
transformation y Q — > y a on the null metric, where = (r, y A ,u), we need 
to know r \, r ^ and r„. From Eq. (53) it follows 



r,x = lv A *VAB,y (54) 



Similarly, 



\u - det(q ABJ C 

'■<■ t y " 

where 



det (UB),C 8 _ 6 

cfei(^) l + ^+p 2y 



^is.c — ( x !ac + x !a x [bc) 9i i + x !a x ^ x % 9i ^ k ( 56 ) 
with the x l £g given functions of (q,p). From Eqs. (53) and (46) 

r,* = in A6 iiAB,* ( 57 ) 

where 

+ 0(A). (58) 



^AB,u 



'I 

A X ,B- 



X a X ~ fi'ij,* 



The null metric r/ a/3 in Bondi coordinates is defined on the world-tube 
Tby 

Note that the metric of the sphere is unchanged by this coordinate trans- 
formation, i.e. r\ AB = fj AB , so we need to compute only the elements rf r , 
if A and rf u on T, or equivalently the metric functions (3, U A and V. From 
Eq. (50), 
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rf r 


= J> 


r 3 f i 


rf A 


= r A 


f A 


rf u 


= r ,& 


fjau . 



5/3 



(r A ) 2 r? AA + 2 r A (r ^ r/ AA - r, fi ) + r ~ A r ^ f] AB 



The contravariant Bondi metric can be written in the form 



-e- 2 ?U 2 

_ e -2/3 





-e-2/3 


r- 2 h 22 r~ 2 h 23 





r -2 ^32 r -2 ^33 












(60) 



(61) 



where Kab is a metric on the sphere of surface area 4tt, such that hABh BC = 
8^ and det{hAB) = det(qAB) = Q, f° r Qab a un h sphere metric. 



4.6. BONDI VARIABLES FOR STARTING UP THE NULL CODE AT THE 
WORLD-TUBE 

The natural variables of the null formalism are certain combinations of 
the null metric functions, which we will give in this section. They will be 
expressed as an expansion in A, as was done in the previous sections, to 
enable us to give these variables in a small region off the world-tube, i.e. at 
points of the grid used to discretize the null equations. These expressions 
are the main results of the extraction module. 

4.6.1. The metric of the sphere J 
Given r and r \, and noting that 



Vab = Vab = r 2 h A B, 



(62) 



h-AB = 
tlAB,\ = 



>VAB 



2r 



^2 I VAB.X 



VAB 



(63) 



In terms of q A , q A and the metric on the unit sphere Jiab, we define the 
metric functions 

J=\q A q B h AB , K=^q A q B h AB , K 2 = 1 + J J. (64) 
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It suffices to evaluate J, since the last relation holds for a Bondi metric. 
We give expressions for J and J \ 

3 = 2^ qAqBvAB 
Jx = ^AV,A-2^ (65) 

Then, in the neighborhood of T, the metric of the sphere is explicitly 
given in Bondi coordinates to second-order accuracy by these last two ex- 
pressions as 

J(y a ) = J + J X X + 0(X 2 ). (66) 



4.6.2. The "expansion factor" (3 

From the last of Eq. (60) we obtain the metric function j3 

/3 = -Ilog(r A ). (67) 

We want to know also its A-derivative, but instead of calculating 7]^ directly 
(which would involve r : \\) 

~-.ru 

/3a = -tt^ = -^ (68) 
we obtain (3^ from the characteristic equation 

P,r = I (j,rJ,r ~ \rf) (69) 

At constant angles (q,p), the relation d\ = r t \d r holds, and we know from 
Eq. (65) J and J a for each outgoing radial null geodesic through the world- 
tube, thus we can write 

A a = (JxJx ~ (Kx) 2 ) (70) 
and from Eq. (64), it follows that 

K X = ±K(JJ X ). (71) 

p^^-^-^mjwr). (72, 

Then, (3 is found to second-order accuracy by: 

/% a ) = /3 + /3 A A + 0(A 2 ). (73) 
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4.6.3. The "shift" U 

The metric function U is related to the Bondi metric, Eq. (61) by 



U = U A q A 



-QA = -\ rT + 



-v 



(74) 



where we have made use of Eq. (60). We also want the A derivative of U 



~XA 



V,x + 



r ,XB r ,B r M \ AB 



r ,x ' 



= " (V + ^ + r f x vf) U + 2 P,x (U + t A n) , (75) 

where in the last line we have used Eq. (68) to eliminate r \\. 
Then, U is found to second-order accuracy by: 



U(y a ) = U + U X X + 0(X 2 ). 



(76) 



4.6.4. The "mass aspect" W 

The metric function V is given in terms of rf r and rf u by 

V = - rr f r /rf u . 



(77) 



For Minkowski-space, V = r at infinity. Thus, we introduce the auxiliary 
metric function W = (V — r)/r 2 , which is regular at null infinity. In terms 
of the contravariant null metric (with the afnne parameter A as the radial 
coordinate) 



W 



1 [rf 



r y T\ 

The A derivative of W is given by 



T ^ T ,B ~AB 

r ,x 



(78) 



W : 



r ,A r ,B ~AB 



fl AD -1 



= - T -fw + ^- [r^ xx + 2 (r A r, XA -r, u ) 
= - r -?((^ + ^)^-K-l)+-J^-rM 



rr a 



(2r A ^ + 2/3, A r A + r^f ) 



r ,A r ,B ~AB 
J2 1 



(79) 
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Then, W is found to second-order accuracy by: 

W(y a ) = W + W : \\ + 0(A 2 ). (80) 

5. Injection 

We have now reached the halfway point in the problem of relating data 
between the Null and Cauchy computational approaches. We turn to in- 
jection, the final part of the problem. Injection is the inverse problem to 
extraction, and allows us to use data obtained by characteristic evolution in 
the exterior of the world-tube in order to provide the exact boundary condi- 
tions for Cauchy evolution in its interior. That is, given the metric in Bondi 
coordinates in a region around the world-tube T, two steps must be carried 
out. The first ingredient is to relate the Cartesian coordinate system to the 
Bondi null coordinates. Secondly, the Cartesian metric components must 
be calculated in the neighborhood of V in the usual way, from the Bondi 
metric components and the Jacobian of the coordinate transformation. 

5.1. LOCATING THE CARTESIAN GRID WORLDLINES 

Given a worldline of a Cartesian grid point, we need to locate, in null 
coordinates, its intersection with a given null hypersurface. Mathematically, 
this means: Given x l and u we want to find t and the null coordinates r, 
as well as y A = y A {x^' 1 ) = (q,p) (where, at this stage, x^ 1 is unknown). 
The starting point is the geodesic equation Eq. (27) 

x a = x i0)a + X£ (0)a + 0(X 2 ), (81) 

Here £^ a = £ a (u, y A (x^ 1 ) is the value of the null vector given by Eq. (26) 
at the point on the world-tube with null coordinates u = f(°), A = and 

We set L a = £ a (u,y A (x i )), so that 

£?(xW)) = L a + 0(\), (82) 

since y A (x l ) = y A (x^ 1 ) + 0(A). (L a can now be interpolated onto the 
stereographic coordinate patch to fourth order accuracy in terms of known 
quantities near the world tube.) 

To the same approximation as (81), we have 

x a = x^ a + \L a + 0(A 2 ) (83) 
so that we have four equations for the five unknowns (x^ l ,t, and A): 

t = n + AL* + 0(A 2 ) (84) 
x i = s (o)i + \V + 0(A 2 ). (85) 
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To find the fifth unknown, we need one additional equation introducing 
new information. This is conveniently supplied by the equation defining 
the world-tube Eq. (20). Even though we have not yet found the values 
for x^' 1 , we may use Eq. (20) to eliminate these unknowns by taking the 
Euclidean norm of Eq. (85), whence using Eq. (19) we find: 

R 2 = f 2 + \ 2 {L ■ L) - 2A(L • x) + 0(A 2 ). (86) 

where the value of R 2 = (x^ ■ x^) is known from the beginning by the 
definition of the world tube, and f 2 = (x ■ x). 

This is a quadratic in A, but we are consistently only keeping terms to 
linear order. This leads to 

A= (f 2 - R 2 )/[2(L-x)] (87) 

Consequently, A = A + 0(A 2 ) and can be calculated in terms of known 
quantities. This leaves four unknown quantities remaining. To proceed to 
evaluate these in terms of known data, we set T = u + 
(X(°),y(°),Z(°)) = x i - KL\ and Y A = (Q,P) = y A (X^ i ). Then, us- 
ing Eq. (18) 

y(°) 

q = rTW)> p = ± rTWv (88) 

on the north (+) and south (— ) patches. T, A and Y A are accurate to 
0(A 2 ) since t = T + 0(A 2 ), A = A + 0(A 2 ) and y A = Y A + 0(A 2 ). 
Finally, the value of r is obtained from 

r = r (°> + rfA, . (89) 

and the values of r a, and r n (used in computing the ff v metric) from 

r,x = rf+r%A, r, A = rf? + r% A, r, u = rg> + rg A (90) 

all to 0(A 2 ) accuracy. 

We calculate r^°\ the derivatives r^, and r$ and the mixed deriva- 
tives r^ A , r^ A and at the points (Q,P) on the world-tube by fourth 
order interpolation. 

The null metric variables are then obtained by fourth order interpolation 
at the Bondi coordinate points (r,Q,P,u), corresponding to the Cartesian 
coordinate points (x l ,T) 

5.2. REVERSING THE COORDINATE TRANSFORMATION 

The following describes the procedure to obtain the contravariant null met- 
ric fj a P in affine coordinates y a = (u, q,p, A) on points on the characteristic 
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grid. This is 90% of the work in obtaining the Cartesian metric at those 
points. What remains is to contract this null metric with the Jacobian 
of the coordinate transformation = x^{y a ). The Jacobian is given in 
Eq. (32) as a series expansion on the parameter A. We will use this ex- 
pression, and consequently we need to work with the contravariant (rather 
than the covariant) metric components, since they will transform with the 
already known Jacobian as: 

9 dy a dy? V ' 1 ' 



5.2.1. The metric on the sphere 

Recall that on each stereographic patch we introduce complex null vectors 
on the sphere q A = (P/2)(5 A + iS A ) where P = 1 + and the q A satisfy 
the orthogonality condition q~AQ A = 2. The q A define a metric on the sphere 
QAB = (4/P 2 )Sab and q A = {2/P){5\ + i5\ ) = q AB q B . 

By construction, the null metric restricted to the sphere, tjab, is the 
same in Bondi coordinates y a and affine coordinates y a . It can be expanded 
as: 

VAb = VAB = cuqAQB + aq A qB + $QAQB + &QAQB (92) 
The values of a and 5 follow from the orthogonality condition 

Aa = VAB q A q B , 45 = r, A Bq A q B = A5, (93) 

From the definition of J and K 

J = ^q A q B h A B, K = ^q A q B h AB , (94) 
and rjAB = r 2 hAB, it follows that 

a = ±r 2 J, 5=^r 2 K, (95) 
so we can write the metric on the sphere in terms of J and K as 

VAB = ^r 2 (J (JaQb + J QAQb + Kq A q B + Kq A q B ) (96) 

From the expressions for the q A , the components of the sphere metric 
are given by: 

v qq = ^r(K + M(j)), n PP = ^-(K-M(j)), Vqp = ^(j) (97) 
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Notice that given the Bondi coordinates of a point, the determinant is 
known 

det(r ?AB ) = r 4 det(q AB ) = \Qr A /P\ (98) 

therefore the inverse sphere metric r/ AB such that rj AB r)BC = &c nas con- 
fidents 









T} qq 


= fj qq = 












= ffP = 








p2 




= f) q P = 


-r^( J )- 



5.2.2. The radial- angular metric coefficients 

We can write the metric coefficients in terms of a single coefficient 7 

fj XA = jq A + jq A (100) 

In components 

fj x P = PM(j), fi Xq = P%(j) (101) 
The value of 7 follows from the expression for the metric function U 

U = U A q A = ^QA = ~ (v XA + r f x V A ^ U, (102) 

2 7 = -U- r -^fj AB q A (103) 
r,\ 

Recall that the derivatives r a, r u are at constant A and r x can be read off 
directly from the Bondi metric function (3 

r x = e" 2/3 (104) 
From the expressions for fj AB and q^, 

r § fj AB q A = |[r 2 (rf p + iff q )+r fl {ff> q + if, m )) 

= ^2^2(K-J) + ir, 3 (K + J)} (105) 
thus the complex coefficient sought is given by 

1 = -\u- e 2 ^ [r 2 (K - J) + ir )3 + J)] (106) 
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The r a can be obtained at the point where the coordinate transformation 
is being performed to 0(A 2 ), by using the angular derivatives of r : \. To this 
end, note that for fixed (q,p): 

r,A = r { °2 + r M \ (107) 

where the derivatives on the right hand side are computed on the world- 
tube. 

5.2.3. The metric function 

By reversing the procedure to obtain W, i.e. from 



W = I [r^ XX + 2 (r A f, xA - r u ) + f, M - 1 j (108) 

we obtain the remaining component of the null metric 

~\\ = e 2(i ^ rW + i_ 2 (r A fj XA - r )U ) - e 2f3 r A r fi fj AB ^j (109) 

We have made use of the relation r \ = e~ 2lS . The r jM and r A are at constant 
A, the former can be evaluated to O(A) accuracy by using the value at the 
world-tube. To get r n to 0(A 2 ), we need to calculate r^ u = — 2e -2 ^/3 u on 

the world-tube, and use r n = r$ + r : \ u \. 

This completes the calculation of the final contravariant metric compo- 
nent, and so the injection process may be carried out at last, after the final 
multiplication by the Jacobian described in Eq. (91). 

5.2.4. Status of the algorithm 

The CCM algorithm described in the previous sections is quite new. It has 
not yet been translated into a stable numerical scheme for black hole col- 
lisions and interactions in 3-D. However, there has already been extensive 
testing of these techniques in simpler situations. Recently, significant codes 
have been built which implement these ideas for scalar waves without grav- 
ity but with full 3-D Cartesian and null spherical coordinate grid patches 
in flat space [17, 18]. Codes have also been written for 1-D problems with 
gravitation, both for spherical self-gravitating scalar waves collapsing and 
forming black holes [19, 20] and for cylindrically symmetric vacuum geom- 
etry [21, 22]. What has not yet been accomplished is the development of 
stable codes with both 3-D geometry and non-trivial strong gravitation, 
although codes currently under development demonstrate the stable evolu- 
tion of a single moving black hole. Thus, we may expect exciting progress 
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in the near future, and the mysterious details of collision between black 
holes may soon be unveiled. 
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